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MULTIDISCIPLINARY 
AERODYNAMIC-STRUCTURAL 
SHAPE OPTIMIZATION USING 
DEFORMATION (MASSOUD) 

Jamshid A. Sa march* 

AVI SA Langley Research Center, Hampton, VA 23681 

This paper presents a multidisciplinary shape parameterization approach. The ap- 
proach consists of two basic concepts: (1) parameterizing the shape perturbations rather 
than the geometry itself and (2) performing the shape deformation by means of the soft 
object animation algorithms used in computer graphics. Because the formulation pre- 
sented in this paper is independent of grid topology, we can treat computational fluid 
dynamics and finite element grids in the same manner. The proposed approach is simple, 
compact, and efficient. Also, the analytical sensitivity derivatives are easily computed 
for use in a gradient-based optimization. This algorithm is suitable for low-fidelity (e.g., 
linear aerodynamics and equivalent laminated plate structures) and high-fidelity (e.g., 
nonlinear computational fluid dynamics and detailed finite element modeling) analysis 
tools. This paper contains the implementation details of parameterizing for planform, 
twist, dihedral, thickness, camber, and free-form surface. Results are presented for a 
multidisciplinary application consisting of nonlinear computational fluid dynamics, de- 
tailed computational structural mechanics, and a simple performance module. 


Nomenclature 

A wing area 

AR wing aspect ratio 

B Bernstein polynomial 

b wing span 

C chord 

Cn I C L coefficients of drag and lift 
c camber 

d degree 

e scale factor for twist and shearing 

N B-spline basis function 

n normal vector 

0 origin of parallelepiped 

P coordinates of NURBS control point 

Q coordinates of the solid elements 

R coordinates of deformed model 

r coordinates of baseline model 

S shearing vector 

T twist plane 

t thickness 

u parameter coordinate 

v MASSOUD design variable vector 

W NURBS weights 

w design variable vector 

X, Y , Z Cartesian coordinates of deformed model 
x, y , z Cartesian coordinates of baseline model 
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a angle of attack, deg 

A total deformation 

deformation 

8 twist angle, deg 

A leading edge sweep angle, deg 

A wing taper ratio 

C f), C coordinates of deformation object 
p twist radius 


Subscripts 

ca camber 

I , J, K total numbers of control points 

i,j,k indices for NURBS control point 

id,jd design variable indices 

in inner 

L wing lower surface 

le leading edge 

m midline 

out outer 

p degree of B-spline basis function in i 

pi planform 

q degree of B-spline basis function in j 

r root 

sh shear 

te trailing edge 

t tip 

th thickness 

tw twist 

U wing upper surface 
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direction 

direction 
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Fuel tank 



Superscripts 

T transpose of the matrix 

Introduction 

M ultidisciplinary design optimization 

(MDO) methodology seeks to exploit the 
synergism of mutually interacting phenomena to 
create improved designs. An MDO process commonly 
involves sizing, topology, and shape design variables. 
Multidisciplinary shape optimization (MSO) finds the 
optimum shape for a given structural layout. Perform- 
ing MSO for a complete airplane configuration is a 
challenging task with high-fidelity analysis tools. The 
analysis models, also referred to as grids or meshes, 
are based on some or all of the airplane components, 
such as skin, ribs, spars, and the stiffeners. The 
aerodynamic analysis uses the detailed definition 
of the skin, also referred to as the outer mold line 
(OML), whereas the computational structural me- 
chanics (CSM) models use all components. Generally, 
the structural model requires a relatively coarse 
grid, but it must handle very complex internal and 
external geometries. In contrast, the computational 
fluid dynamics (CFD) grid is a very fine one, but 
it only needs to model the external geometry. The 
MSO of an airplane must treat not only the wing 
skin, fuselage, flaps, nacelles, and pylons, but also 
the internal structural elements such as spars and 
ribs (see Fig. 1). The treatment of internal structural 
elements is especially important for detailed finite 
element (FE) analysis. For a high-fidelity MSO 
process to be successful, the process must be based 
on a compact and effective set of design variables 
that yields a feasible and enhanced configuration. 
For more details, readers are referred to an overview 
paper by this author on geometry modeling and grid 
generation for design and optimization. 1 

Model parameterization is the first step for an 
MSO process. Over the past several decades, shape 
optimization has been successfully applied for two- 
dimensional and simple three-dimensional configura- 
tions. Recent advances in computer hardware and 
software have made MSO applications more feasible 
for complex configurations. An important ingredient 
of aerodynamic shape optimization is the availability 
of a model parameterized with respect to the aerody- 


namic design variables, such as planform, twist, shear, 
camber, and thickness. The parameterization tech- 
niques can be divided into the following categories: 
discrete, polynomial and spline, computer-aided de- 
sign (CAD), analytical, and deformation. Readers are 
referred to reports by Haftka and Grandhi, 2 Ding, 3 
and Samareh 4 for surveys of shape optimization and 
parameterization. 

In a multidisciplinary application, the parameter- 
ization must be compatible with and adaptable to 
various analysis tools ranging from low-fidelity tools, 
such as linear aerodynamics and equivalent laminated 
plate structures, to high-fidelity tools, such as nonlin- 
ear CFD and detailed CSM codes. Creation of CFD 
and CSM grids is time-consuming and costly for a full 
airplane model: detailed CSM and CFD grids based 
on a CAD model take several months to develop. To 
fit the MSO process into product development cycle 
times, the MSO must rely on parameterizing the anal- 
ysis grids. For a multidisciplinary problem, the process 
must also use a geometry model and parameterization 
consistently across all disciplines. Gradient-based op- 
timization requires accurate sensitivity derivatives of 
the analysis model with respect to design variables. 

This paper presents a shape parameterization ap- 
proach suitable for MSO as part of a multidisciplinary 
design optimization application. The approach con- 
sists of two basic concepts. The first concept is based 
on parameterizing the shape perturbations rather than 
the geometry itself. The second concept is based on 
using the soft object animation 5 (SOA) algorithms 
for shape parameterization. The combined algorithm 
initially introduced by this author 1 was successfully 
implemented for aerodynamic shape optimization with 
analytical sensitivity for structured grid 6,7 and un- 
structured grid 8 CFD codes. This algorithm has also 
been used for multidisciplinary application of a high- 
speed civil transport (HSCT). 9 ’ 10 

Parameterizing the Shape Perturbations 

At first sight, parameterization by splines may seem 
to be a viable approach for shape parameterization. 
The spline representation uses a set of control points, 
which could define any arbitrary shape. These control 
points could be used as design variables for optimiza- 
tion. Typically, over a hundred control points are 
required to define an airfoil section, and over 20 airfoil 
sections to define a conventional wing. This require- 
ment results in over two thousand control points (i.e., 
six thousand shape design variables) for a simple wing. 
The number of control points is even larger for a com- 
plete airplane model created with a commercial CAD 
system. The large number of control points is needed 
more for accuracy than for complexity. 

Even if using a large number of design variables were 
feasible, the automatic regeneration of analysis models 
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(e.g., CSM and CFD grids) is not possible with cur- 
rent technology. For example, it takes several months 
to create an accurate CSM model of an airplane. Also, 
traditional shape parameterization processes parame- 
terize only the OML and are ineffective in parameteriz- 
ing internal components such as spars, ribs, stiffeners, 
and fuel tanks (see Fig. 1). 

It is possible to use any shape (e.g., a sphere) as 
the initial wing definition, allowing the optimizer to 
find the optimum wing shape; however, this option is 
not commonly practiced. Typically, the optimization 
starts with an existing wing design, and the goal is 
to improve or redesign the wing performance by us- 
ing numerical optimization. The geometry changes 
(perturbations) between initial and optimized wing are 
very small, 11,12 but, the difference in wing performance 
can be substantial. An effective way to reduce the 
number of shape design variables is to parameterize 
the shape perturbations instead of parameterizing the 
shape itself. Throughout the optimization cycles, the 
analysis grid can be updated as 


R(v) — r + AR(v) (1) 

For multidisciplinary aerodynamic-structural shape 
optimization using deformation (MASSOUD), the 
change, A R, is a combination of changes in thickness, 
camber, twist, shear, and planform: 


A R — SRth + 3Rca + 3Rtw + 3Rsh + 3Rpi (2) 

Far fewer design variables are required to parameterize 
the shape perturbations A R than the baseline shape 
r itself. 

Figures 2 and 3 depicts the typical MSO and MAS- 
SOUD processes. In a typical MSO process (Fig. 2), 
a geometry modeler perturbs the baseline geometry 
model. Because automatic grid generation tools are 
not available for all disciplines, automating this MSO 
process would be very difficult. In contrast, the MAS- 
SOUD process (Fig. 3) relies on parameterizing the 
baseline grids and avoids the grid generation process, 
hence automating the entire MSO process. 

Soft Object Animation 

The held of SOA in computer graphics 5 provides 
algorithms for morphing images 13 and deforming mod- 
els. 14,15 These algorithms are powerful tools for mod- 
ifying the shapes: they use deformation of a high-level 
shape, as opposed to manipulation of lower level ge- 
ometric entities. Hall presents an algorithm and pro- 
vides computer codes for morphing images. 13 The de- 
formation algorithms are suitable for deforming mod- 
els represented by either a set of polygons or a set of 
parametric curves and surfaces. The SOA algorithms 
treat the model as rubber that can be twisted, bent, 



Fig. 2 A typical MSO process. 



Fig. 3 The MASSOUD process. 


tapered, compressed, or expanded, but will still retain 
its topology. This technique is ideal for parameteriz- 
ing airplane models that have external skin as well as 
internal components (e.g., see Fig. 1). The SOA al- 
gorithms link vertices of an analysis model (grid) to a 
small number of design variables. Consequently, the 
SOA algorithms can serve as the basis for an efficient 
shape parameterization technique. 

Barr presented a deformation approach in the con- 
text of physically based modeling. 14 This approach 
uses physical simulation to obtain realistic shape and 
motions and is based on operations such as transla- 
tion, rotation, and scaling. With Barr’s algorithm, 
the deformation is achieved by moving the vertices of 
a polygon model or the control points of a paramet- 
ric curve and surface. Sederberg and Parry presented 
a variant, 15 of the free form deformation (FFD) algo- 
rithm, which operates on the whole space, regardless of 
the representation of the deformed objects embedded 
in the space. The algorithm allows a user to rnanip- 
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ulate the control points of trivariate Bezier volumes. 
The disadvantage of FFD is that the design variables 
may have no physical significance for the design en- 
gineers. This drawback makes it difficult to select an 
effective and compact set of design variables. This pa- 
per presents a set of modifications to the original SOA 
algorithms to alleviate this and other drawbacks. 

For the modified SOA algorithms presented in the 
next several sections, implementation will include the 
following common set of steps: 

1. Select an appropriate deformation technique and 
object. This step defines the forward mapping 
from the deformation object coordinate system 
(£; . T j. ( j to the baseline grid coordinate system 
(x,y,z). 

2. Establish a backward mapping from the baseline 
grid coordinate system ( x , y, z ) to the deforma- 
tion object coordinate system (£, i), ( ) . The £, i), ( 
mapping parameters are fixed and are indepen- 
dent of the shape perturbations. This preprocess- 
ing step is required only once. 

3. Perturb the control parameters (design variables) 
of the deformation object. 

4. Evaluate the grid perturbation (A R) and shape 
sensitivity derivatives (dR/dv) with the £, rj, pa- 
rameters. 

The following sections provide recipes for using SOA 
algorithms for parameterizing airplane models for 
thickness, camber, twist, shear, and planform changes. 

Thickness and Camber 

We used a nonuniform rational B-spline (NURBS) 
representation as the deformation object for thickness 
and camber parameterization. The NURBS represen- 
tation combined the desirable properties of National 
Advisory Committee for Aeronautics (NACA) defini- 
tion 16 and spline techniques, and it did not deteriorate 
nor destroy the smoothness of the initial geometry. 

The changes in thickness and camber were repre- 
sented by 

>': V;,(i)>: .Y.,( V )U; ; 

V) = — — j- ^7 (3) 

y~l Ni tP (£,) y~] 

i= 0 j= 0 

>: V;, (i)>: v,,( f/ )iu ; 
a a (e, V) = ^7 (4) 

i = 0 j = 0 

Figures 4 and 5 show the NURBS control points in 
(£,??) and (x, y, z) coordinate systems, respectively. 


% 



Fig. 4 Thickness and camber definitions in wing 
coordinate system. 



Fig. 5 Thickness and camber definitions in x, y, 
and ^ coordinate system. 


The control points and weights could be used as design 
variables. 

The NURBS representation had several important 
properties for design and optimization. A NURBS 
curve of order p. having no multiple interior knots, 
is p — 2 differentiable. As a result, the NURBS rep- 
resentation was able to handle a complex deformation 
and still maintain smooth surface curvature. Read- 
ers are referred to the textbook by Farm for details on 
the properties of NURBS representation. 1 ‘ The control 
points were the coefficients of the basis functions, but 
the smoothness was controlled by the basis functions, 
not by the control points. The NURBS representa- 
tion was local in nature, allowing the surface to be 
deformed locally, hence leaving the rest of the surface 
unchanged. Equations (3-4) served as the forward 
mapping between the thickness and camber design 
variables and the grid perturbation (SRth,3Rca)- 

The next step was to establish the backward map- 
ping from the deformation object (i.e., the NURBS 
surface) coordinates (£, rf) to the baseline model coor- 
dinates ( x,y,z ). The percentage of chord, %C, was 
used for £, and the spanwise location, y, was used for 
V- 

£ = %C, n = y (5) 

To calculate %C , we needed to determine the wing 
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the B-spline basis functions. 



TO p 


Fig. 6 Curves defining the backward mapping. 


chord at each y station. The baseline CAD model 
provided the leading edge Ri e { y), trailing edge Rt e (y), 
wing midline R m (y ) , and normal vector defining the 
airfoil plane T(i )), as shown in Fig. 6. The curve defin- 
ing the wing midline did not have to be at the center 
of the wing, but needed to be somewhere between the 
upper and the lower wing surfaces. The Ri e (y), Rt e [y), 
and R m (y) were used to separate points on the upper 
surface from points on the lower surface. 

Because we knew y for each grid point, we were able 
to define a plane that passed through the grid point 
by means of a normal vector defined by T(y). We then 
used the following equations to find the intersection of 
this plane and the curves shown in Fig. 6. 


f(y)-[r-R le (y)] T = 0 

(6) 

f(^)-[r- J R te (ij)] T = 0 

(7) 

T(y) ■ [r- R m (v)] T = 0 

(8) 


Equations (6-8) must be solved for each grid point in 
the model. For a high-order NURBS curve, Eqs. (6- 
8) are nonlinear and can be solved by the Newton- 
Raphson method. The solution to Eqs. (6-8) for each 
y was a set of three points located at the leading edge, 
the trailing edge, and the center. The %C was cal- 
culated based on the leading and trailing edge points. 
Next, we needed to separate the grid points defining 
the wing model into upper and lower surfaces. We 
connected the three points obtained from Eqs. (6- 
8) to form a curve that separated the upper surface 
from the lower surface. This curve need not represent 
the camber line accurately, and a wing with drooping 
leading edge or with highly cambered airfoil sections 
may require more than one R m \y) to define the curve. 
With this approach, the deformation may be localized 
to a specific design area by setting allowable %C m i n , 

Umax : ^minj and 7] max - 

As the design variables (control points Pij) 
changed, we calculated the contribution to the thick- 
ness and camber by Eqs. (3-4). The advantage of this 
process was that the sensitivity of grid point location 
with respect to design variables was only a function of 


dR _ dR 
dPth i d:jd dP caid 


Njd,p (QNjd,q(ll)WidJd 

E N iA€)T, N j,q( r i)Wij 

i—0 j = 0 

(9) 


Consequently the sensitivity, as shown in Eq. (9), was 
independent of the design variables ( Pid,jd ) and the 
coordinates ( x , y, z). Thus, we needed to calculate the 
sensitivity with respect to thickness and camber only 
at the beginning of the optimization. 


Twist and Shear 

The twist angle is defined as the difference between 
the airfoil section incident angle at the root and each 
airfoil section incident angle. Similarly, the shear (di- 
hedral) is defined as the difference between the airfoil 
leading edge z coordinate for the root and the z co- 
ordinate at each airfoil section. If the twist angle at 
the tip is less than the twist at the root, the wing is 
said to have a washout, which could delay the stall at 
the wing tip. Also, as the wing washout increases, the 
wing load shifts from outboard to inboard. As a re- 
sult, the spanwise distribution of the twist angle plays 
an important role in the wing performance. 

The SOA algorithms are used to modify the wing 
twist and shear distribution. Alan Barr presented a 
series of SOA algorithms for twisting, bending, and 
tapering an object. 14 Watt and Watt referred to these 
algorithms as nonlinear global deformation. 5 Seder- 
berg and Greenwood extended Barr’s ideas to handle 
complex shapes. 18 Modified versions of these algo- 
rithms are presented in this paper. 

To modify the twist and shear distributions, the 
wing was embedded in a nonlinear deformation ob- 
ject referred to as a twist cylinder, shown in Fig. 7. 
The center of the cylinder was defined by a NURBS 
curve, R m (y). The effect of deformation was confined 
to a section of a wing by limiting the parameter y to 
vary between 7^ m ; n and ?? max . The y m - m could be ex- 
tended to the wing root, and the ^ max went beyond the 
wing tip. The cylinder could be twisted and sheared 
only in a plane (twist plane) defined by a point along 
R m (y) with a normal vector of T(rj). The pi n {v) and 
Pout{r)) were the radii of the inner and outer cylinders, 
respectively (see Fig. 7). The deformation had no ef- 
fect on grid points located outside the outer cylinder, 
and the effect of deformation was scaled linearly from 
the outer cylinder to the inner cylinder. This linear 
blending allowed us to blend the deformed region with 
the undeformed region in a continuous manner. 

The 8(y) and S(y) variables are defined by the 
NURBS representations: 
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Fig. 7 Inner twist cylinder. 


J2 N i,p{v)Wi6i 

^ ( 10 ) 

E N i>p (r,)Wi 

i = 0 

E N iiP (T,)WiSi 

—J ( 11 ) 

E NiAriW 

i = 0 

Similarly to thickness and camber algorithms, we used 

V = y, T{ri) — (0, y, 0) T (12) 

The second step for twist and shear deformation was 
to establish the forward mapping from the deformation 
object (twist cylinder) coordinate system ( r ) ) to the 
model coordinate system ( x,y,z ). We used Eq. (8) to 
determine r). Once y was determined, we calculated 
the local p(rj), p m {i]), p ou t{rj), T(r)), 8{rj), and S{rj). 
The point r was rotated 8(i)) deg about R m ,(i)) and 
sheared S(i)). 

SR tvl (rj) = e(p)p(r))[sm8(r)),0,cos8(i])] T (13) 

6R*h{v) = e(p)S(r ) ) (14) 

where e(i)) was a scale factor that diminished the effect 
of deformation as we approached the outer cylinder. 



Fig. 8 Twist definition for a transport. 



tip. 



Fig. 10 Planform of a generic HSCT. 

cos 6(i)) depended on the twist design variables and 
were updated for every cycle of the optimization. In 
contrast, the term dS(i))/dSi was independent of shear 
design variables Si (see Eq. (11)). 

Figure 8 shows the inner twist cylinder for a com- 
mercial transport. Figure 9 shows the result of twist- 
ing the wing 45 deg at the tip. This amount of twist is 
large and unrealistic, but demonstrates the effective- 
ness of the SOA. 


%) = 
Sfa) = 


{ 0 if p{v) > Pont iv) 

ll-7cE if Pm - < p™ t(i) ( 15 ) 

1 if p{r}) < Pin 

The sensitivity of a grid point with respect to the 
twist and shear design variables was 


dR 

Wi 

dR 


5(p)/o(T)^^[ cos0 ( , 7) ; °.- sin6 '( , ?)]I 16 ) 


5 (p) 


dS(rj) 

dSi 


(17) 


Planform Parameterization 

The wing planform is typically modeled with a set of 
two-dimensional trapezoids in the x-y plane. Figure 10 
shows the planform of a generic HSCT that uses two 
trapezoids. As shown in Fig. 11, each trapezoid is 
defined by the root chord ( C r ), tip chord ( Ct ), span 
(6), and sweep angle (A). From these values, other 
planform parameters such as area (A), aspect ratio 
(AR), and taper ratio (A), are defined: 

A=|(C r + a), A = g- (18) 


The term d6(i))/d9i was independent of the twist de- The FFD algorithm described by Sederberg and 

sign variables 9{ (see Eq. (10). However, sin #(77) and Parry 15 is ideal for deforming the polygonal models. 
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Fig. 11 Planform definition. 




Like other SOA algorithms, this algorithm maintains 
the polygon connectivity, and the deformation is ap- 
plied only to the vertices of the model. The FFD 
process is similar to embedding the grid inside a block 
of clear, flexible plastic (deformation object) so that, 
as the plastic is deformed, the grid is deformed as well. 
Deformation of complex shapes may require several 
deformation objects. The shapes of these deformation 
objects are not arbitrary. In fact, the shapes must be 
three-dimensional parametric volumes and could range 
from a parallelepiped as shown in Fig. 12 to a general 
NURBS volume as shown in Fig. 13. The block is 
deformed by perturbing the vertices that control the 
shape of the deformation block (e.g., corners of the 
parallelepiped). For parametric volume blocks, param- 
eters controlling the deformation are related through 
the mapping coordinates (£,??, C)- These coordinates 
are used in both forward and backward mapping. 


Figure 12 shows a general parallelepiped defined by 
a set of control points forming three primary edges 
or directions along £, rj. and £. The relation for a 
parallelepiped is defined as 


r(i,riX) = 0 + + n v r) + n ( ( (19) 

where O is the origin of the parallelepiped, and ft f, n Tj . 
and n.f are the unit vectors along the parallelepiped 
primary edges in the £, r ) , and ( directions, respec- 
tively. Equation (19) defines a mapping between 
the deformation object (parallelepiped) and the grid 
points. The grid points, r, are mapped to the coordi- 
nates of the parallelepiped, £, ??, and ( , as 


n-r, x ng • (r - P 0 ) 
n , , x • (fif) 
n.£ x h( ■ (r - Tp) 
fj£ x • (n v ) 
n( x n v ■ (r - Pp) 
Jl( x hri ■ (fif) 


( 20 ) 


A grid point is inside the parallelepiped if 0 < £, r), ( < 

1 . 

The FFD technique based on the parallelepiped is 
very efficient and easy to implement. This technique 
is suitable for local and global deformation. The only 
disadvantage is that the use of the parallelepiped lim- 
its the topology of deformation. To alleviate this 
disadvantage, Sederberg and Parry proposed to use 
nonparallelepiped objects. 15 They also noted that the 
inverse mapping would be nonlinear and would require 
significant computations. 

Another popular method to define FFD is to use 
trivariate parametric volumes. Sederberg and Parry 
used a Bezier volume. 15 Coquillart at INRIA ex- 
tended the Bezier parallelepiped to nonparallelepiped 
cubic Bezier volume. 19 This idea has been fur- 
ther generalized to NURBS volume by Lamousin and 
Waggenspack. 20 The NURBS blocks are defined as 


E Nj, P 2{n) E N kiP3 {QWi tj ,kPi,j,k 

i = 0 j- 0 fc = 0 

E^,pi(0 E N jip2 ( V ) t N ktP3 (QW iJ>k 

i = 0 ?' = 0 k — 0 

( 21 ) 

Lamousin and Waggenspack 20 used multiple blocks to 
model complex shapes. This technique has been used 
for design and optimization by Yeh and Vance 21 and 
also by Perry and Balling. 22 

The common solid elements used in FE analysis 
(Fig. 14) can be used as deformation objects as well. 


r(Z,Tl, C) 
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Fig. 14 FE analysis solid elements. 


The mapping from the solid element coordinates is de- 
fined 23 by 

= ( 22 ) 

i 

where N,- are the FE basis functions and Qi are the 
nodal coordinates of deformation objects, which are 
related to the design variables. The equations for in- 
verse mapping are nonlinear for all solid elements with 
the exception of tetrahedron solid elements. The solid 
elements provide a flexible environment in which to de- 
form any shape. Complex shapes may require the use 
of several solid elements to cover the entire domain. 

To model the planform shape, we used hexahedron 
solid elements with four opposing edges parallel to the 
z coordinate. Then, the planform design variables 
were linked to the corners of the hexahedral elements. 
Figure 15 shows the initial and deformed models for 
a transport configuration. The solid lines represent 
the controlling hexahedron solid elements. The base- 
line model is on the left-hand side, and the deformed 
shape is on the right-hand side. 

As with the camber and thickness algorithms, the 
sensitivity of grid point coordinates was independent 
of the design variables (P) and coordinates ( x,y,z ). 
Thus, we needed to calculate the sensitivity only once, 
at the beginning of the optimization. 

Implementation 

Figure 16 shows the implementation diagram for the 
combined algorithm. The implementation started with 
a CAD model that defined the geometry. The first two 
steps were implemented in parallel. The first step was 
to determine the number and the locations of the de- 
sign variables with the aid of the CAD model. In the 
second step, the grids were manually generated for all 
involved disciplines. In the third step, the forward and 
backward mappings described in the previous sections 
were calculated for each grid point. In the fourth step, 
the new grid was deformed in response to the new 
design variables, and the sensitivity derivatives were 



Fig. 15 Planform deformation of a transport. 



Fig. 16 MASSOUD Implementation diagram. 


computed as well. The third and fourth steps were 
completely automated. The first three steps were con- 
sidered preprocessing steps and needed to be done only 
once. 

Parameterizing CSM Models 

Parameterizing CFD and CSM models appears to 
be similar in nature, but the CSM model parameter- 
ization has two additional requirements. First, the 
CSM model parameterization must include not only 
the OML but also the internal structural elements, 
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such as spars and ribs. Second, the deformed CSM 
model must be a valid design. For example, the spars 
must stay straight during the optimization. The algo- 
rithms presented in this paper can easily handle the 
first requirement. However, if the planform design 
variables are not selected with care, the second re- 
quirement could easily be violated. To avoid creating 
an invalid CSM model, the model must be parame- 
terized with few hexahedron solid elements, and those 
used must be aligned with major structural compo- 
nents such as spars and ribs. 

Sensitivity Analysis 

Sensitivity derivatives are defined as the derivatives 
of the coordinate locations with respect to the design 
variables. The previous sections present a formula- 
tion for shape parameterization based on a specific 
set of design variables (c; , i. = l,* max ). It is pos- 
sible to introduce a new set of design variables (wj, 
j = 1 , jrmix) - The sensitivity derivatives with respect 
to Wj were computed based on the chain rule differen- 
tiation as 


OR OR dvi 

dwj dvi dwj 


(23) 


The previous sections provide techniques to compute 
the first term on the right-hand side. The second term 
is defined in a matrix form where the matrix has i max 


rows and j max columns. 



r dv-\ 
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dv i "I 


dwi 

dw 2 

dw 7rT1 


dv-2 

dv 2 

dv 2 


dw 1 

dw -2 



8 »>m ax 

dv tmax 
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dwi 

dw 2 
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Design-Variable Sequencing 

In a typical optimization problem, the number of 
design variables is determined a priori. However, the 
NURBS representation algorithms permits the use of 
an adaptive algorithm to determine the number of de- 
sign variables. The design variables are control points 
of a NURBS curve or surface. Optimization of a wing 
section could start with as few as three design variables 
(see Fig. 17). Then, the number of design variables can 
be increased to five by enriching the NURBS curve or 
surface. Enrichment is accomplished by inserting ad- 
ditional knots into the design NURBS curve or surface. 

This method is similar to mesh sequencing and 
multigrid methods used in CFD to accelerate the con- 
vergence. Multigrid methods exhibit a convergence 
rate independent of the number of unknowns in the 
discretized system. 24 Using MASSOUD for design- 
variable sequencing needs to be investigated further. 



3 DVs 


5 DVs 


9 DVs 


17 DVs 


33 DVs 


Fig. 17 A sequence of design- variable sets. 



Fig. 18 Baseline and deformed CSM models of an 
HSCT. 


Results and Conclusions 

The algorithms presented in this paper have been 
applied for parameterizing a simple wing, a blended 
wing body, and several HSCT configurations. Fig- 
ure 18 shows the baseline and deformed CSM grids for 
an HSCT. The solid lines represent the initial position 
of the hexahedron solid elements controlling the plan- 
form variation. The parameterization results from this 
research have been successfully implemented for aero- 
dynamic shape optimization with analytical sensitivity 
for structured 6 and unstructured 8 CFD grids. 

An aerodynamic optimization of an ONERA M6 
wing was performed 6 using a sequential linear pro- 
gramming technique. The objective of the optimiza- 
tion was to minimize the drag while maintaining the 
same lift as the baseline design. Figure 19 shows the 
design cycle history for both lift and drag. In this 
optimization, the angle of attack is fixed, and it was 
found that in order to move away from the current 
design, the constraint on the lift coefficient had to be 
relaxed temporarily. This is shown clearly in the fig- 
ure: for the first 19 design cycles, Cl was allowed to 
deviate by up to 0.01 from the desired value. After 
design cycle 19 the tolerance on the lift constraint was 
tightened to 10 -6 . The net result was approximately 
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Fig. 19 Design cycle history of ONERA M6 wing 
optimization for coefficient of drag. 

29 counts of drag reduction at the baseline lift. Fig- 
ure 20 shows comparisons of the solutions computed 
on the baseline and final designs. The results indicate 
a significant reduction in the shock strength at most 
spanwise stations. 

This approach has also been applied to multidisci- 
plinary optimization of a HSCT. 9 ’ 10 Figure 21 shows 
the design cycle history for aircraft drag, as measured 
relative to the baseline values. The figure shows the 
drag has been reduced by 7.5 % relative to the base- 
line. Although the optimizer has not fully converged 
for this case, the convergence history from 20 design 
cycles suggests that little additional reduction in drag 
would be obtained from additional design cycles. Fig- 
ure 21 shows a comparison of the baseline and final 
surface pressures on both the upper and lower sur- 
faces of a HSCT. The planform changes that occurred 
between the initial and final design cycles are also rep- 

AIAA- 



Fig. 20 Comparison of the M6 wing final design 
and baseline surface pressures. 



Design Cycle 

Fig. 21 Design cycle history of HSCT optimization 
for drag. 

resented. The primary effect on the planform has been 
to increase the span and aspect ratio slightly and to 
move the outer wing leading edge break to a more in- 
board spanwise location. Although not evident in the 
figure, the wing thickness has been slightly reduced. 

The parameterization algorithm presented in this 
paper is easy to implement for an MDO application 
for a complex configuration. The resulting parame- 
terization is consistent across all disciplines. Because 
the formulation is based on the SOA algorithms, the 
analytical sensitivity is also readily computed. The 
algorithms are based on parameterizing the shape 
perturbations, thus enabling the parameterization of 
complex existing analysis models (grids). Another 
benefit of parameterizing the shape perturbation is 
that the process requires few design variables. Use of 
NURBS representation provides strong local control, 
and smoothness can easily be controlled. 
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Fig. 22 Comparison of the HSCT final design and 
baseline surface pressures. 
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